/*
DO FILE FOR REPLICATION FOR 
Local labor market effects of nuclear power plants

BY Duha T. Altindag, Reem El Cheikh Taha, Jennifer U. Jones, and R. Alan Seals, Jr., May 2025
*/

**"${path}clean_data/Final Panel.dta" includes all counties in the USA
**"${path}clean_data/Considered panel.dta" includes only the counties that constructed their first reactor between 1974 and 1978 which is our sample

clear all
set more off

global path "C:/Users/Reem/Desktop/stata files/"

********************FIGURES

*******Figure 1: The Number of Nuclear Reactor Construction Licenses Awarded Over Time
***CL
use "${path}clean_data/power stations.dta", clear
bys CP_ISSUED_year: egen num=count(UNIT) if CP_ISSUED_year!=.
collapse (mean) num,by(CP_ISSUED_year)
ren CP_ISSUED_year year
drop if year==.
save "${path}temp_data/number of reactors getting CL per year.dta", replace

*******COL
clear
use "${path}clean_data/power stations.dta", clear
bys COMMERCIAL_OPERATION_year: egen numb=count(UNIT) if COMMERCIAL_OPERATION_year!=.
collapse (mean) numb,by(COMMERCIAL_OPERATION_year)
ren COMMERCIAL_OPERATION_year year
drop if year==.
save "${path}temp_data/number of reactors getting COL per year.dta", replace

*******Merge together
use "${path}temp_data/number of reactors getting COL per year.dta", clear
merge 1:1 year using "${path}temp_data/number of reactors getting CL per year.dta"
drop _merge
replace num=0 if num==.
replace numb=0 if numb==.
sort year
set obs 66
// Define the list of years and the corresponding observation numbers
local years 1991 1992 1994 1995 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2014 2019 2020
local obs   46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
// Loop through each pair and apply the replacement
forvalues i = 1/21 {
    local y : word `i' of `years'
    local o : word `i' of `obs'
    replace year = `y' in `o'
}
replace num=0 if num==.
replace numb=0 if numb==.
***Graph
gen top=numb+num
	twoway ///
    (bar top year, barwidth(0.65) color(maroon)) ///
    (bar num year, barwidth(0.65) color(navy)) ///
    , ///
    ytitle("The Number of Licenses") xtitle("") ///
    xlabel(1960(10)2020) ///
    legend(order(2 "Construction License" 1 "Commercial Operation License"))
*manually change the intensity of the maroon and navy bars to 100% instead of the default 80%
graph export "${path}results/Figure1.png", replace



********Figure 2: Time Between the Year of Announcement and Issuance of Commercial License
use "${path}clean_data/power stations.dta", clear
drop if ACTION=="future" |ACTION=="cancelled"
gen start= ANNOUNCED_year
gen end= COMMERCIAL_OPERATION_year
gen duration=.
replace duration=end-start
egen id=group(CP_ISSUED_year)
sort id
bys id: egen av_duration=mean(duration)

preserve
keep id CP_ISSUED_year av_duration
duplicates drop
twoway ( connected av_duration CP_ISSUED_year, sort), xtitle("Year When Construction Started") ytitle("Duration Between Announcement and Operations (Years)") xscale(range(1950 1980)) xlabel(1950(5)1980)
graph export "${path}results/Figure2.png", replace
restore

*******Figure 3: Map
**states and county coordinates
use "${path}clean_data/Final Panel.dta",clear
keep if county_considered==1
keep GEOFIPS treated  
duplicates drop
merge 1:1 GEOFIPS using "${path}clean_data/cb_2019_us_county_20m/us_data.dta"
spmap treated using "${path}clean_data/cb_2019_us_county_20m/us_coords.dta" if STATEFP!="02" & STATEFP!="15" & STATEFP<="56", id(geoid) ///
    fcolor(Greys) legend(off)
***fix plot 1 to Grayscale15 in the graph editor
***fix plots 2 and 3 same border color
graph export "${path}results/Figure3.png", replace

*********Figures 4 & 5: ES
use "${path}clean_data/Considered Panel.dta", clear

egen g=group(GEOFIPS  year)
bys g:gen plant_exists=construction_dummy+commercial_dummy
replace plant_exists=1 if plant_exists>=2

***Emp
gen treat_year= year if plant_exists==1
bys GEOFIPS: egen first_treat=min(treat_year)
gen ry=year-first_treat
gen never_treat= (first_treat==.)
g before5= ry<=-5 //5 years and before
forval j=4(-1)1 {
g before`j'= ry==-`j'
}
forval j=0/19 {
g post`j'= ry==`j'
}
g post20= ry>=20 //20 years and after

areg emp_rate before5-before2 post0-post20    i.year, a(GEOFIPS) cluster(GEOFIPS)
preserve
clear
set obs 26
g t=_n-6
drop if t==-1
g b=.
g se=.
forval j=2/5 {
replace b=_b[before`j'] if t==-`j'
replace se=_se[before`j'] if t==-`j'
}
forval j=0/20 {
replace b=_b[post`j'] if t==`j'
replace se=_se[post`j'] if t==`j'
}
g ub=b+1.96*se
g lb=b-1.96*se
label define label_t -5 "-5+" 20 "20+"
label values t label_t
twoway ///
    (rcap lb ub t if t < 0, lcolor(maroon)) /// 
    (rcap lb ub t if t >= 0, lcolor(maroon)) /// 
    (scatter  b t if t < 0, mcolor(navy) msymbol(circle)) /// 
    (scatter  b t if t >= 0, mcolor(navy) msymbol(circle)) /// 
    ,xlabel(-5 0 5 10 15 20, labels valuelabel) xtitle("Year Relative to the Construction Start") ytitle("Estimate") ///
      legend(off)
	  graph export "${path}results/Figure4.png", replace
restore

***Wages
areg log_wages_Rpercapita before5-before2 post0-post20  i.year, a(GEOFIPS) cluster(GEOFIPS)
preserve
clear
set obs 26
g t=_n-6
drop if t==-1
g b=.
g se=.
forval j=2/5 {
replace b=_b[before`j'] if t==-`j'
replace se=_se[before`j'] if t==-`j'
}
forval j=0/20 {
replace b=_b[post`j'] if t==`j'
replace se=_se[post`j'] if t==`j'
}
g ub=b+1.96*se
g lb=b-1.96*se
label define label_t -5 "-5+" 20 "20+"
label values t label_t
twoway ///
    (rcap lb ub t if t < 0, lcolor(maroon)) /// 
    (rcap lb ub t if t >= 0, lcolor(maroon)) /// 
    (scatter  b t if t < 0, mcolor(navy) msymbol(circle)) /// 
    (scatter  b t if t >= 0, mcolor(navy) msymbol(circle)) /// 
    ,xlabel(-5 0 5 10 15 20, labels valuelabel) xtitle("Year Relative to the Construction Start") ytitle("Estimate") ///
      legend(off)
	  graph export "${path}results/Figure5.png", replace
restore

******Figures 6 & 7 Sun & Abraham 
use "${path}clean_data/Considered Panel.dta", clear
egen g=group(GEOFIPS  year)
bys g:gen plant_exists=construction_dummy+commercial_dummy
replace plant_exists=1 if plant_exists>=2

***Emp
preserve
gen treat_year= year if plant_exists==1
bys GEOFIPS: egen first_treat=min(treat_year)
gen ry=year-first_treat
gen never_treat= (first_treat==.)
g before5= ry<=-5 //5 years and before
forval j=4(-1)1 {
g before`j'= ry==-`j'
}
forval j=0/19 {

g post`j'= ry==`j'
}
g post20= ry>=20 //20 years and after

eventstudyinteract emp_rate before5-before2 post0 post1-post20, cohort(first_treat) control_cohort(never_treat) absorb(i.GEOFIPS i.year)
matrix C = e(b_iw)
mata st_matrix("A",sqrt(st_matrix("e(V_iw)")))
matrix C = C \ A
matrix list C
coefplot matrix(C[1]), se(C[2]) yline(0, lcol(red)) vertical xtitle("Year Relative to the Construction Start") ytitle("Estimate") xlabel(1 "-5+" 5 "0" 10 "5" 15 "10" 20 "15" 25 "20+" , labels valuelabel) 
graph export "${path}results/Figure6.png", replace
restore

***Wages
preserve
gen treat_year= year if plant_exists==1
bys GEOFIPS: egen first_treat=min(treat_year)
gen ry=year-first_treat
gen never_treat= (first_treat==.)
g before5= ry<=-5 //5 years and before
forval j=4(-1)1 {
g before`j'= ry==-`j'
}
forval j=0/19 {

g post`j'= ry==`j'
}
g post20= ry>=20 //20 years and after

eventstudyinteract log_wages_Rpercapita before5-before2 post0 post1-post20, cohort(first_treat) control_cohort(never_treat) absorb(i.GEOFIPS i.year)
matrix C = e(b_iw)
mata st_matrix("A",sqrt(st_matrix("e(V_iw)")))
matrix C = C \ A
matrix list C
coefplot matrix(C[1]), se(C[2]) yline(0, lcol(red)) vertical xtitle("Year Relative to the Construction Start") ytitle("Estimate") xlabel(1 "-5+" 5 "0" 10 "5" 15 "10" 20 "15" 25 "20+" , labels valuelabel) 
graph export "${path}results/Figure7.png", replace
restore


******Figures 8: ES Employment Industry Specific
use "${path}clean_data/Considered Panel.dta", clear
egen g=group(GEOFIPS  year)
bys g:gen plant_exists=construction_dummy+commercial_dummy
replace plant_exists=1 if plant_exists>=2
gen treat_year= year if plant_exists==1
bys GEOFIPS: egen first_treat=min(treat_year)
gen ry=year-first_treat
gen never_treat= (first_treat==.)
g before5= ry<=-5 //5 years and before
forval j=4(-1)1 {
g before`j'= ry==-`j'
}
forval j=0/19 {

g post`j'= ry==`j'
}
g post20= ry>=20 //20 years and after

***construction industry
areg const_emprate before5-before2 post0-post20 i.year, a(GEOFIPS) cluster(GEOFIPS)
preserve
clear
set obs 26
g t=_n-6
drop if t==-1
g b=.
g se=.
forval j=2/5 {
replace b=_b[before`j'] if t==-`j'
replace se=_se[before`j'] if t==-`j'
}
forval j=0/20 {
replace b=_b[post`j'] if t==`j'
replace se=_se[post`j'] if t==`j'
}
g ub=b+1.96*se
g lb=b-1.96*se
label define label_t -5 "-5+" 20 "20+"
label values t label_t
twoway ///
    (rcap lb ub t if t < 0, lcolor(maroon)) /// 
    (rcap lb ub t if t >= 0, lcolor(maroon)) /// 
    (scatter  b t if t < 0, mcolor(navy) msymbol(circle)) /// 
    (scatter  b t if t >= 0, mcolor(navy) msymbol(circle)) /// 
    ,xlabel(-5 0 5 10 15 20, labels valuelabel) xtitle("Year Relative to the Construction Start") ytitle("Estimate") ///
      legend(off)
graph save "${path}results/const_emp.gph", replace
restore

***utilities industry
areg trans_emprate before5-before2 post0-post20 i.year, a(GEOFIPS) cluster(GEOFIPS)
preserve
clear
set obs 26
g t=_n-6
drop if t==-1
g b=.
g se=.
forval j=2/5 {
replace b=_b[before`j'] if t==-`j'
replace se=_se[before`j'] if t==-`j'
}
forval j=0/20 {
replace b=_b[post`j'] if t==`j'
replace se=_se[post`j'] if t==`j'
}
g ub=b+1.96*se
g lb=b-1.96*se
label define label_t -5 "-5+" 20 "20+"
label values t label_t
twoway ///
    (rcap lb ub t if t < 0, lcolor(maroon)) /// 
    (rcap lb ub t if t >= 0, lcolor(maroon)) /// 
    (scatter  b t if t < 0, mcolor(navy) msymbol(circle)) /// 
    (scatter  b t if t >= 0, mcolor(navy) msymbol(circle)) /// 
    ,xlabel(-5 0 5 10 15 20, labels valuelabel) xtitle("Year Relative to the Construction Start") ytitle("Estimate") ///
      legend(off)
graph save "${path}results/utilities_emp.gph", replace
restore

***gov industry
areg gov_emprate before5-before2 post0-post20 i.year, a(GEOFIPS) cluster(GEOFIPS)
preserve
clear
set obs 26
g t=_n-6
drop if t==-1
g b=.
g se=.
forval j=2/5 {
replace b=_b[before`j'] if t==-`j'
replace se=_se[before`j'] if t==-`j'
}
forval j=0/20 {
replace b=_b[post`j'] if t==`j'
replace se=_se[post`j'] if t==`j'
}
g ub=b+1.96*se
g lb=b-1.96*se
label define label_t -5 "-5+" 20 "20+"
label values t label_t
twoway ///
    (rcap lb ub t if t < 0, lcolor(maroon)) /// 
    (rcap lb ub t if t >= 0, lcolor(maroon)) /// 
    (scatter  b t if t < 0, mcolor(navy) msymbol(circle)) /// 
    (scatter  b t if t >= 0, mcolor(navy) msymbol(circle)) /// 
    ,xlabel(-5 0 5 10 15 20, labels valuelabel) xtitle("Year Relative to the Construction Start") ytitle("Estimate") ///
      legend(off)
graph save "${path}results/gov_emp.gph", replace
restore

graph combine "${path}results/const_emp.gph" "${path}results/utilities_emp.gph" "${path}results/gov_emp.gph"
graph export "${path}results/Figure8.png", replace

*******Figure 9: ES Wages Industry Specific
use "${path}clean_data/Considered Panel.dta", clear
egen g=group(GEOFIPS  year)
bys g:gen plant_exists=construction_dummy+commercial_dummy
replace plant_exists=1 if plant_exists>=2
gen treat_year= year if plant_exists==1
bys GEOFIPS: egen first_treat=min(treat_year)
gen ry=year-first_treat
gen never_treat= (first_treat==.)
g before5= ry<=-5 //5 years and before
forval j=4(-1)1 {
g before`j'= ry==-`j'
}
forval j=0/19 {

g post`j'= ry==`j'
}
g post20= ry>=20 //20 years and after

***construction industry
areg const_earning_Rpercapita before5-before2 post0-post20 i.year, a(GEOFIPS) cluster(GEOFIPS)
preserve
clear
set obs 26
g t=_n-6
drop if t==-1
g b=.
g se=.
forval j=2/5 {
replace b=_b[before`j'] if t==-`j'
replace se=_se[before`j'] if t==-`j'
}
forval j=0/20 {
replace b=_b[post`j'] if t==`j'
replace se=_se[post`j'] if t==`j'
}
g ub=b+1.96*se
g lb=b-1.96*se
label define label_t -5 "-5+" 20 "20+"
label values t label_t
twoway ///
    (rcap lb ub t if t < 0, lcolor(maroon)) /// 
    (rcap lb ub t if t >= 0, lcolor(maroon)) /// 
    (scatter  b t if t < 0, mcolor(navy) msymbol(circle)) /// 
    (scatter  b t if t >= 0, mcolor(navy) msymbol(circle)) /// 
    ,xlabel(-5 0 5 10 15 20, labels valuelabel) xtitle("Year Relative to the Construction Start") ytitle("Estimate") ///
      legend(off)
graph save "${path}results/const_earning_Rpercapita.gph", replace
restore

***utilities industry
areg trans_earning_Rpercapita before5-before2 post0-post20 i.year, a(GEOFIPS) cluster(GEOFIPS)
preserve
clear
set obs 26
g t=_n-6
drop if t==-1
g b=.
g se=.
forval j=2/5 {
replace b=_b[before`j'] if t==-`j'
replace se=_se[before`j'] if t==-`j'
}
forval j=0/20 {
replace b=_b[post`j'] if t==`j'
replace se=_se[post`j'] if t==`j'
}
g ub=b+1.96*se
g lb=b-1.96*se
label define label_t -5 "-5+" 20 "20+"
label values t label_t
twoway ///
    (rcap lb ub t if t < 0, lcolor(maroon)) /// 
    (rcap lb ub t if t >= 0, lcolor(maroon)) /// 
    (scatter  b t if t < 0, mcolor(navy) msymbol(circle)) /// 
    (scatter  b t if t >= 0, mcolor(navy) msymbol(circle)) /// 
    ,xlabel(-5 0 5 10 15 20, labels valuelabel) xtitle("Year Relative to the Construction Start") ytitle("Estimate") ///
      legend(off)
graph save "${path}results/utilities_earning_Rpercapita.gph", replace
restore

***gov industry
areg gov_earning_Rpercapita before5-before2 post0-post20 i.year, a(GEOFIPS) cluster(GEOFIPS)
preserve
clear
set obs 26
g t=_n-6
drop if t==-1
g b=.
g se=.
forval j=2/5 {
replace b=_b[before`j'] if t==-`j'
replace se=_se[before`j'] if t==-`j'
}
forval j=0/20 {
replace b=_b[post`j'] if t==`j'
replace se=_se[post`j'] if t==`j'
}
g ub=b+1.96*se
g lb=b-1.96*se
label define label_t -5 "-5+" 20 "20+"
label values t label_t
twoway ///
    (rcap lb ub t if t < 0, lcolor(maroon)) /// 
    (rcap lb ub t if t >= 0, lcolor(maroon)) /// 
    (scatter  b t if t < 0, mcolor(navy) msymbol(circle)) /// 
    (scatter  b t if t >= 0, mcolor(navy) msymbol(circle)) /// 
    ,xlabel(-5 0 5 10 15 20, labels valuelabel) xtitle("Year Relative to the Construction Start") ytitle("Estimate") ///
      legend(off)
graph save "${path}results/gov_earning_Rpercapita.gph", replace
restore

graph combine "${path}results/const_earning_Rpercapita.gph" "${path}results/utilities_earning_Rpercapita.gph" "${path}results/gov_earning_Rpercapita.gph"

graph export "${path}results/Figure9.png", replace







	

